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1. Introduction 

Two separate themes from the Computer Vision literature come together in this paper: the use of na- 
tionally symmetric operators, and the idea that several modules of visual perception require that the "most 
conservative 1 ' solution that meets a given set of boundary conditions be computed. The two themes are 
combined in an investigation of which operator to use in the interpolation of smooth surfaces from one- 
dimensional boundary constraints. Such constraints arise naturally in a variety of visual problems. 

In the next section we review the role of rotationally symmetric operators in Computer Vision, and we 
derive conditions which linear and quadratic forms in the first and second directional derivatives must satisfy 
in order to be rotationally symmetric. We then discuss the idea that vision is a conservative process, citing 
examples from both figure perception and scene analysis. The "most conservative" solution is modelled using 
the calculus of variations to find the minimum function that satisfies a given performance index. A major 
problem associated with the use of the calculus of variations is guaranteeing the existence of an minimum 
fimction (see for example Courant and Hilbert 1953, p.173). A theorem of Grimson(1981, theorem 2) proves 
that a sufficient condition for the existence of a minimum is that the performance index should be a seminorm 
on the space of functions. The condition is not necessary. For example, Horn(1981) has determined the curve 
that minimizes die integral square curvature subject to tangency conditions at the end points; the performance 
index is not a seminorm. 

Grimson(1981) notes that many intuitively plausible performance indices based on mean and Gaussian 
curvature are not seminorms, but that the square Laplacian f 2 xx + 2f xx f yy + f 2 yy and the quadratic variation 
fix + y 2 xy + fly are - We show here that <my quadratic form in f XX9 f xyi and f yy is a seminorm. 

To further constrain the choice of performance index in the infinite set of quadratic forms, we require 
in addition that the quadratic form should be rotationally symmetric. We prove that there are essentially two 
different choices: the square Laplacian and the quadratic variation. All the remaining possibilities are linear 
combinations, that is, form a vector space with these two as a basis. 

To choose between the square Laplacian and the quadratic variation, we consider their respective Ruler 



conditions and natural boundary conditions (Courant and Hilbert, 1953). The Euler conditions are identical, 
but the natural boundary conditions, which are derived from the statics of a deformed thin plate, favor the 
quadratic variation since they offer tighter constraint in this case. 

2. Rotationally symmetric operators in vision 

A major concern of Computer Vision is the isolation of constraints that combine with the information 
provided in the image to yield an interpretation. Early work on polyhedra (Clowes 1971, Huffman 1971, 
Mackworth 1973, Waltz 1972, Sugihara 1978, 1981, Kanade 1981) focussed upon the discovery of constraints 
deriving from the image forming process, constraints that relate image fragments, like junctions and lines, to 
their scene counterparts, vertices and edges. As Computer Vision turned its attention away from plane-faced 
objects to the natural world, other constraints were required. Often the constraints expressed some facet of the 
intuitive notion of "smoothness" and did so in a way that supported useful computations (Strat 1979, Brooks 
1979, Ikeuchi and Horn 1981, Woodham 1978, Horn and Schunck 1981). Recently, smoothness and image 
forming have been combined using differential geometry (Grimson 1981, Witkin 1981, Binford 1981). 

One constraint that is usually implicit, but is occasionally made explicit, expresses the idea that perceptual 
processes are often' approximately isotropic. It seems that humans usually do not show strong directional 
preferences when detecting edges, motion, or reflectance boundaries. We seem to be equally adept at per- 
ceiving the layout and orientation of a visible surface regardless of its orientation relative to the view vector. 
Ullman(1976) argues for an explicit isotropy constraint in his work on subjective contours (see also Knuth 
1979). 

Processes that are isotropic are naturally computed by rotationally symmetric operators, since the values 
they return are unaffected by the coordinate system chosen for the image. Conversely, rotationally symmetric 
operators compute isotropic information. As we shall see, many operators that have been proposed for vision 
are not rotationally symmetric but directionally selective. Some authors have, however, proposed rotationally 
symmetric operators, particularly for early visual processing. 
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Precise definitions of rotational symmetry for functions, operators (or functional), and, by specialization, 
matrices are given in the following section. In the rest of this section we assume that the definitions are already 
understood. 

Some kinds of blurring in an image forming system can be approximated by convolution with a 
Gaussian. The rotationally symmetric Gaussian can be defined by: 

1 — r 2 

G(r)==~7ra 2 exp(-^). 

Pratt(1978) presents several techniques, such as convolution with the generalized inverse of the blur 
function, for restoring the image, (see for example, his figures 14.2.1, 14.3.2). 

The Laplacian A = f xx + f yy is well known to be rotationally symmetric* and its use has been proposed 
several times in Computer Vision and Image Processing. If an image is blurred in a way that can be ap- 
proximately modelled by passing the image through a system with a Gaussian point spread function, then it 
can be sharpened by subtracting a multiple of its Laplacian (Rosenfeld and Kak 1976, p.184, Prewitt 1970, p. 
107). Pratt(1978, figure 17.4.5) illustrates the use of the Laplacian for enhancing the edges in an image. 

Weska, Dyer and Rosenfeld(1976) note that convolving a step edge with a Laplacian operator gives rise to 
a pulse pair: a negative pulse at the transition from the lower plateau to the edge, and a positive pulse at the 
transition from the edge to the upper plateau (see also Horn 1974, Marr and Hildreth 1980). They suggested 
that the image intensities at the locations of the positive and negative pulses could be used to set thresholds to 
use in segmenting the image into regions. 

Several authors have noted the relative insensitivity of human perception to small intensity gradients 
(Hcrskovits and Binford 1970, Marr 1976, Marr and Hildreth 1980, McCann et al. 1974). They have noted 
that the effect can be explained by assuming that the vision system uses operators approximating second 
derivatives. This so-called lateral inhibition effect seems to be performed by center surround operators in 
the retina (see for example Richter and Ullman 1980). The Laplacian is a rotationally symmetric second 

t A proof of this is given in Section 3 below. 



differential operator, and an attractive candidate to perform lateral inhibition. 

The use of the Laplacian for edge detection was proposed by Horn(1974) in a study of the determination 
of lightness. Following Land and McCann(1971), Horn restricted attention to images of planes colored with 
patches of uniform reflectance or color. Within a patch, grey level variations are due to small variations 
in illumination, and they are smooth compared to the abrupt changes between patches. The conventional 
approach to detecting significant changes in intensity had been to note that the gradient of the image is 
small within a region, but is infinite across a reflectance boundary between regions. For a particular image 
tesselation and quantization of grey levels, the gradient is always finite. It is usually much larger, however, at 
a reflectance boundary than it is within a region. Horn(1974) rejected using the gradient since "the first partial 
derivatives are directional and thus unsuitable since they will for example completely eliminate evidence of 
edges running in a direction parallel to their direction of differentiation." The Laplacian is the lowest order 
linear combination of derivatives that is rotationally symmetric. A reflectance boundary can be detectedby the 
paired positive and negative peaks on either side of the boundary, and localized by noting the position where 
the Laplacian crosses zero between the peaks^. 

Marr and Hildreth(1980) have proposed that edges are detected in the human visual system by an 
operator that approximates AG, where A is the Laplacian, and G is a rotationally symmetric Gaussian. We 
shall show in the next section that the application of a rotationally symmetric operator, such as the Laplacian, 
to a rotationally symmetric function, such as the Gaussian, is itself rotationally symmetric. It follows that the 
Marr-Hildreth operator is rotationally symmetric. Marr and Hildreth note that intensity changes occur at a 
number of scales and are often superimposed. They suggest that an image should be smoothed by a number 
of bandpass filters to isolate the changes at a particular range of scales. The Gaussian is chosen as the filter to 
optimize localization of changes in both the spatial and frequency domains. 

We noted above that the Gaussian and the Laplacian have figured prominently in early visual processing. 
The Gaussian has mostly been used to approximate the point spread function corresponding to the blurring of 



t See Binford(1981) for more on the distinction between detection and localization of an intensity change. 
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a point source. Marr and Hildreth deliberately introduce Gaussian blurring. They further note that AG can be 
approximated by a difference of Gaussians, G\ — G2. Nishihara and Larson(1981) note that the difference of 
Gaussians is to be preferred on grounds of efficiency. Macleod(1972) proposes an edge detection operator that 
is the difference of two Gaussians. However, no analysis of its performance is given, and no indication is given 
that the operator approximated a low-pass filtered second derivative. 

Regarding the use of the Laplacian, Marr and Hildreth do not seem to make isotropy an explicit con- 
straint on edge detection. Instead, Hildreth(1980,page 13) notes that "a number of practical considerations, 
which will be illuminated in the discussion of the implementation, suggested that the . . . operators not be 
directional". Suppose instead that directional operators are used. The simplest algorithm for edge detection 
has two stages. First, the image is convolved with the directional operators in "sufficiently many" directions. 
Second, the outputs are combined to determine the orientation and extent of intensity changes. Regarding 
the first stage, both Marr and Hildreth(1980, page 193) and Hildreth(1980, page 40) claim that the cost of 
convolving the image with a "sufficient" number of operators is excessive. They show that a single rotationally 
symmetric operator (the Laplacian) gives precisely the same results if a condition called "linear variation" 
holds. Regarding the second stage, Hildreth(1980, page 36) observes that edges in a direction close to that of 
the mask are elongated in the direction of the mask. She also notes that operators at several orientations give 
significant responses to any given edge, and that combining the responses is non-trivial. 

There are two essentially different issues here that need to be clearly separated. Intensity changes first 
have to be detected and then localized as a set of "feature points" marking the position of the change in the 
image, and characteristics of the corresponding edge. The detection of feature points is inherently isotropic, 
as Horn(1974) noted. The feature points have then to be combined to produce descriptions of edge segments. 
Edge segments are clearly directional, indeed a central problem concerns the determination of the direction of 
an edge in an image. The computation of rich descriptions of edge segments is, as Hildreth notes, not at all 
easy. Marr's(1976) original Primal Sketch work was almost entirely concerned with it. Binford( 1981) discusses 
the application of directional operators to compute the directionality of an edge. 



The Gaussian and Laplacian are not the only rotationally symmetric operators that have been proposed 
in computer vision. Prewitt(1970, p. 107) observes that "derivatives of all orders can be used to form isotropic 
nonlinear differential operators, provided that derivatives of odd order appear only in even functions. The 
simplest of these ... is the squared gradient", namely V • V, where V is the column vector 
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Earlier in the same article, Prewitt(1970, p. 85) suggests that "the Hankel transformation enters naturally 
in the analysis of systems with isotropic point spread functions and greatly facilitates restoration." The sugges- 
tion does not appear to have been investigated in computer vision. 

We noted earlier that an important aspect of modelling perception is the isolation of constraints which 
capture some facet of smoothness. Horn and Schunck(1981) consider the determination of optical flow fields 
and note that "if every point of the brightness pattern can move independently, there is little hope of recover- 
ing the velocities". One way to express the additional constraint of smoothness is to minimize the integral of 
the performance index 

S(u,ti)-(uH«}) + '(^ + «5), 

where u and v are the x and y components of the optical flow, and subscripts denote partial differentiation. 
We shall show in the next section that this operator is rotationally symmetric. In many simple situations the 
smoothness constraint is significantly wrong only at occluding boundaries. 

We conclude this review of the use of rotationally symmetric operators in vision with Grimspn*$(1981) 
work on surface interpolation. As it will be the focus of Section 5, our remarks will be brief. The Marr-Poggio 
theory of human stereo vision yields the disparity (scaled depth) at matched edge points that are computed 
by the Marr-Hildreth approach described above. The disparity map is as sparse as the set of matched edge 
points, whereas human perception is of smooth surfaces passing through the given disparity points. Crimson 
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(1981) interpolates a smooth surface from the given set of edge points by a local parallel algorithm that applies 
a rotationally symmetric operator to minimize the quadratic variation introduced above. 

3. Conditions for rotational symmetry 

A function fM 2 h-+ 3? is rotationally symmetric if its polar form is only dependent on radial distance 
r = (x 2 + y 2 )^ and not on direction <j> = tan"" 1 1. Clearly, a function is rotationally symmetric if and only 
if it can be represented as a function of (x 2 + y 2 )i. An alternative definition can be given that is often more 
convenient for functions, and that can be generalized to operators. A function is rotationally symmetric if and 
only if it yields the same value under an arbitrary rotation of coordinates. 

An anticlockwise rotation from one set of image coordinates (x } y) to another (X, Y) is effected by a 
rotation matrix: 



cos<£ sin^> 
— sin <f> cos 0. 



= R 



(0) 



For convenience, we shall denote cos<£ by c and sin^ by s. To simplify notation, we shall not make 
explicit the dependence of the rotation matrix R on the angle <f>. A function / is rotationally symmetric if and 
only if the untransformed version f(x, y) is equal to the transformed version f(X, Y). We shall occasionally 
find it useful to borrow the mathematical shorthand that equates a function f(X, Y) with a function of a single 

rp 

vector argument f(R[x,y] ). 

Example 1. The function fi(x, y) = (x 2 -f- y 2 ) is rotationally symmetric: 

f l (X,Y) = ((xc + ys) 2 + (yc-x8) 2 ) 
= (x 2 + y 2 ) 

= M x > y)- 
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Example 2 The function fi{x, y) = xy is not rotationally symmetric: 



h[X,Y) = {xc + ys){yc-xs) 

yl — X 2 

= zycos2<^4- — 5 — sin2^, 

it 

and so f 2 (X, Y) = f 2 {x, y) only when <l> = or <£ = *. 

We can extend the definition of rotational symmetry to operators 

0:(3ft 2 k» 3ft) h* (SR 2 .-+ »). 

An operator is rotationally symmetric if 0(/) is a rotationally symmetric function, for all functions 
fM 2 ^9l. 

Example 3. The function produced by the operator 0^ defined by 

OiC/)(*,y) = e«**) 

is rotationally symmetric if and only if/ is. In general then, the operator Oi is not rotationally symmetric. 
However, the Gaussian is a rotationally symmetric operator as it combines examples 1 and 3. 

Most of the operators of interest in computer vision are combinations of the first and second directional 
derivatives £, ^, gi, ^, djj&, and gp. We need to determine the effect of a coordinate rotation on these 
directional derivatives. By the chain rule, 

dx ~*~ dxdX + SxdY 

_ d d 

~ C dX~ 8 dY- 

Similarly, 



d_ d_ d_ 
dy~ 8 dX +C dY 



It follows that 



= R T 



a* 



where T denotes matrix transpose. Since R is a rotation matrix, its transpose equals its inverse, so 



a*" 



= R 



(1) 



Operators in general, and differential operators in particular, depend upon the choice of coordinate 
frame. To make the dependence of the differential operator on the choice of coordinate frame explicit, we 
introduce the notation 



j^***i. 



7 (*,y)- 



With this notation, equation (1) becomes 



V ( x,y)=i*V, 



(*.«)> 



(2) 



where ^( x ,y) is the column vector 



Proposition 1. Linear combinations of ^ and §q are not rotationally symmetric. 
Proof. Any linear form in the first directional derivatives has the form 

[* /4 v (*,v)- 



ITic condition for rotational symmetry is 
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By equation (2), 

and so the linear differential operator is rotationally symmetric if and only if 

[X fi] = [X ti]R t 

so that [X p] is an eigenvector of R. The eigenvalues of J? are c + is and c — is. So there are no real 
eigenvectors unless </> is a multiple of ?r. Since the condition is not satisfied for all <f>, no linear combination is 
rotationally symmetric, i 

The same style of analysis can be applied to other combinations of first derivatives such as the operator 

It is easy to show that 2 (x,y) is not equal to 2 ( X)y y for example when <j> — f . ■ 
In section 2, we referred to an operator proposed by Prewitt(1970), namely 



UJ'+UJ 1 



that is, the vector dot product 



V(x,y)V(*>y)* 



More generally, we often consider quadratic differential expressions such as 
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V f*,v) 



y i. 



'(*,v)- 



Such an expression is called -a quadratic form if the matrix is symmetric, that is \i = ja By equation 1, 



so that 



ifandonlyif 



V ( x,y)=#V 



(*iV> 



V (^) MV (*,y) = V fx,r) MV (W 



i^Mfl = M, 



where/? is an arbitrary rotation matrix, and 



M = 



b e 



Since the transpose R T of a rotation matrix R is the inverse of R, a quadratic form is rotationally 
symmetric if and only if the corresponding matrix M commutes with all rotation matrices. We will refer to 
matrices M having this property as being rotationally symmetric. 

Lemma 1. A 2 by 2 matrix is rotationally symmetric if and only if it has the form 



M = 



X fi 

— M \ 
Proof. We require RM = MR for all rotation matrices H, that is 
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"c 


— 8 


'X fl 




"X \£ 


*c 


—a" 


8 


c m 


y I 




y t 


8 


c m 



Expanding, and equating terms, this holds if and only if 

X = £. 

Alternatively, only the operations of scaling by a constant k and multiplication by a rotation matrix R f 
commute with all rotation matrices in two dimensions So M — kR! for some scale factor k and some rotation 
matrix R!% 

Proposition 2. Up to scaling, the only rotationally symmetric quadratic form in £ and ^ is V {Xj y) • V (a ., y ). 

Proof. A quadratic form in £ and ^ has the form 



y (*,y) 



X f£ 



V 



(x,V> 



(3) 



To be rotationally symmetric, as well as symmetric (so that it is a quadratic form), Lemma 1 implies that 



H — 0. 

It follows that the matrix in equation (3) is X/ 2 .| 

The operator f 2 x -f f y is commonly used as a measure of the contrast across an intensity change. Notice 
that other popular measures of the contrast, such as (f x -f fy) 2 , (/* — fyf, or \\f x \\ + ||/ y || are not rotationally 
symmetric, and therefore respond differently to edges in different directions (see Rosenfeld and Kak 1976, 
p279). 

We now consider linear and quadratic forms in £?, £%, ggj, and &p. It is convenient to not assume 
^-j = ^j for the developments that follow. - 

TTic first task is to find a matrix R* so that 
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a? 

dzBy 

& 

By&x 



du 1 J 



= #* 



dm? 

. w 1 . 



(4) 



The (i, j) element of the matrix J? t will be denoted by ry. Applying the chain rule as before, but this 
time to relate the second derivatives in (X, Y) to those in (x, y), we find that the four by four matrix R* can be 
written in the form 



R* = 



r n R T r 2l R T 
r n R T r 22 R T 



(5) 



Definition 1. (ben Israel and Greville 1974, page 41)Let/4 = [a y ] andB = [6y] be m by m and n by n 
matrices respectively. The mn by ran matrix A ® S, called the Kronecker product of A andB, is defined by 
multiplying each element a(£, j) of .4 by the matrix B, to form the block matrix 



'a n B 


CL12B . 


. aimP 


a 2 \B 


&22JB 


• ^m^ 


flm\B 


a m 25 * 


Q"mmP 



(6) 



With this notation, 



R*=R T ®R T , 



/**N 



so that 



t Recall the definition of the matrix R from equation (0). 
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i?* = 



C 2 — 8C 
SC C 2 



8C 8 
2 



-8 



8C —8 2 C 2 



—8C 
— 8C 



(7) 



U' 8C SC <T 

Note that the elements of A <g) J3 are naturally indexed by 4-tuples: 



{i4 ® J5> WJbl = 

We state without proof a number of. simple properties of the <g) operation. They are essentially 
straightforward consequences of the properties of ordinary multiplication, and are stated without proof. 

Proposition 3 

■(•) {A<g>B) T = A T <8)B T 

(it) (i4®B)- 1 =i4r 1 ®B-- 1 . 

. m {iii) {A®B)<g)C = A<g){B®C) 

For the remainder of the paper, we restrict attention to the application of to -R and its transpose. 
Proposition 4. The rotationally symmetric linear combinations of J^, g^, ^, and J^ are linear 
combinations oftheLaplacian A = £5 + jjfa, and the smoothness measure ^ — 3^. 
Proof. Let the linear combination be 



[x m v i\ 



& 

Bx&y 
&y5x 

.<h? J 



Following the proof of Proposition 1, the condition for rotational symmetry is 



fX M y. 6]# r <g)H r -[X P M Q, 



f\ 
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for all rotation matrices R and the corresponding rotation angle </>. Expanding R T <g) R T by equation (7), we 
find 



-SC —SC 8* 



so that 



[\ n v i) 



8C C* 



-5* —8C 



8C — 8 2 C 2 



-8C 

LB" 8C 8C C 2 J 



[* f* v g], 



/**s 



It follows that 



[^ i* v i\ 



— s 2 — sc — sc s 2 



SC — 8* — S" 6 — 8C 



sc —s 2 —a 2 — sc 
S 2 sc sc — S 2 J 



[0 0]. 



[i — $ **+*/ o o] 



—2s 2 — 2sc 

2sc -2s 2 





The determinant of the upper left 2 by 2 submatrix is 



— 2sc 


2s 2 " 


-2s 2 


— 2sc 











. 



= [000 o]. 



(4s 4 + 4s 2 c 2 ) = 4a 2 . 

Since this is not zero for all angles <£, X — £ and \i + v arc both zero. A basis for the infinite set of 
linear combinations satisfying these conditions is provided by setting X and \i equal to one, which proves the 
Proposition. | 

Before turning to quadratic forms, analogous to Proposition (2), we define a projection, operator on 

R T (g) R T that makes explicit the assumption f xy = f yx . 
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Definition 2. Let D = \d i3 ] be a 4 by 4 matrix. The projection of D is the 3 by 3 matrix D* : 

rfll (dl2+<*13) dl4 

{<k\ + <fei) (d 2 2 + d 32 + eta + ^33) (^24 + ^34) 

CJ41 (^42 + ^43) ^44 

That is, the second and third columns as well as the second and third rows are combined by addition. 
Proposition 5. 

[abb c]D[ a b b cf 



is equivalent to 

[a b c)D [a b c) t 

where £>* is the projection of D. 

The proof is by equating terms, and is omitted. We now give the main result of this section. 

Proposition 6. The rotationally symmetric quadratic forms in ^, ^^, ^, and ^ form a vector 
space. If ^ = ^, the matrices associated with the rotationally symmetric quadratic forms project to 3 by 
3 matrices of the form 



2a 

P <* + /? 

It follows that the rotationally symmetric quadratic forms that satisfy ^j == ^ form a vector space 
that has the quadratic variation, ( J^) 2 + 2(^) 2 + (|^) 2 , and the square Laplacian, (^ + jg>) 2 , as a 
basis. 

Proof. Since the matrix in a quadratic form is defined to be symmetric, a quadratic form in ^, ^|^, 

*&> and ^ can be written 



/"""V 



^•"N 
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#5 



6x 



f x dv 2 



\A B* 
IJB T C m 


SxSy 
dybx 



where A and C are symmetric 2 by 2 matrices, and B is 2 by 2. As usual, the quadratic form is rotationally 
symmetric if and only if 



R T ®R? 



A B 



\B T C\ \B T C. 



A B 



R T ®R T , 



where R is an arbitrary rotation matrix. It follows that 



cR T sR r 
sR T cR T . 



A B 



A B' 
\B T C. 



cR T sR 7 ! 
l-sR T cR T 



and hence that 



cR T A + sR T B T cR T B + sR T C ' 
- S R T A + cR T B T - 8 R T B + cR T C 



' cAR T - sBR T sAR T + cBR T 
cB T R T -sCR T aB T R T + cCR T 



Equating submatrices, we find that for all rotation angles <f> 



/»**N 



c{R T A - AR T ) + s{R T B T + BR T ) = 0, 

c(R T C - CR T ) - s(B T R T + R T B) = 0, 

s(CR T - R T A) + c(R T B T - B T R T ) = 0, 

s{R T C - AR T ) + c{R T B - BR T ) = 0. 

Consider equation (10) or (11) when <j> = f . Equating terms, we find that 

an = 422 

Ol2 = — 321 
&l\ — — ^12- 



(8) 

(9) 

(10) 

(11) 



(12) 
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Similarly, equation (8) or (9) when $ — § yields 



&ll+&22 = 0- 



(13) 



Expanding equation (8) for general <j> yields 

bu +a 12 *=0, (14) 

622—021 = 0, (15) 

621 + &12 + «22 — an = 0. (16) 

Combining equations (12) through (16) we find that in order to be rotationally symmetric, the matrix 



A B' 
\JB T C 
has the form 



a + /? 7 -7 /* 

7 a — 6 6 7 

— 7 6 a — 6 — 7 

L P 7 ~7 « + # 
A matrix of this form projects to 



~a + P P 
2a 

/? a + /3 

where a = 612 — an and /3 = 612. It is easy to show that linear combinations of matrices of this form are 
of the same form, so that the rotationally symmetric quadratic forms constitute a vector space. Clearly, the 
square Laplacian and the quadratic variation, corresponding to the cases a = 1,/? = and a .— 0,/3 = 1 
respectively, form a basisj 



r^ 
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We show that the measure of smoothness of optical flow proposed by Horn and Schunck(1981) is rota- 
tionally symmetric. Recall from section 2 that the measure is defined by the operator 

S(u,v) = ( U l + uD + (vl + vl). 

We extend the Kronccker product operator <g) to vectors, and then show how to define S(u, v) in terms 
of vector Kronecker products. 

Definition 3. (a) Let a = [a\. . .a m ] and b = [61. . .6 n ] be vectors. The Kronecker product of a and b is the 
ran dimensional vector [a\bi. . .a\b n 0261. . .a m 6 n ], 

(b) By extension, if = [O x . . .O m ] is a vector of operators and/ = [/1. . ./ n ] is a vector of functions, the 
Kronecker product of and / is the ran dimensional vector of functions 

[OiCA)-.^i(/„)...o m c/ n )]. 

The components u and t; of optical flow are functions of x t y, and f. Recall that ^( x ,y) = [£% ^) T - 
According to definition 3, 



v ( ., w) (8)iu «r-i s ^ s ^i, 



so that 



S(u, t;) = (V ( , itf) <g> [u *] T ) • (v ( ,, y) <g> [u t,] T ). 



If the coordinate frame is rotated through <j> by the matrix R, the optical flow components become R[u v] T . 
The Horn-Schunck measure is rotationally symmetric if and only if 

(R®Rf(R®R)=I 4f 
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where J 4 is the 4 by 4 identity matrix. The rotational symmetry is a simple consequence of Proposition 3. 
A rotationally symmetric operator has the general form 

( s, y) (V,V®V,V®V®V,H, 
and its application to a rotationally symmetric function f{x f y) has the form 

{X) y){f{x t y)). 

To see that this is rotationally symmetric, we rotate the coordinate frame to (X, Y) by a matrix R as before. 
Since and / are rotationally symmetric, all the occurences of R (including its Kronecker square, cube, and so 
on) introduced by the frame change can be deleted. It follows that the application of a rotationally symmetric 
operator to a rotationally symmetric function is itself rotationally symmetric. In particular, the A(G) filters of 
the Marr-Hildreth theory of edge detection are rotationally symmetric. 

4. Vision as a conservative process 

The second theme of this paper is that a number of vision modules construct the most conservative inter- 
pretation that is consistent with the given data, and that is subject to an appropriate set of suitably formulated 
constraints. A major concern of Computer Vision has always been the isolation of constraints that enable the 
interpretation of an image. Constraints embody observations about the way the world is, at least, most of the 
time. Although such observations can be as specific as cataloging familiar figures and shapes, it has proved 
more fruitful to first uncover constraints that correspond to general observations that are widely applicable. 
Constraints are used together with the data computed from the image to construct an interpretation. The 
representations of the information from the image and the constraints determine, and are determined by, 
the interpretation process. For example, early blocks world programs represented constraints as catalogs of 
labcllings, an approach that led naturally to search processes for interpretation (Clowes 1971, Kanade 1981). 
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As Computer Vision has turned its attention to images of the natural world, constraints have concerned 
the smoothness of surfaces and movement. The relationship to boundary value problems of physics and 
mathematics suggests itself. The information computed from the image sets the boundary conditions, and the 
constraints determine which (and whether a) solution to the boundary value problem is found. Horn(1974) 
solved an instance of Poisson's problem using Green's functions to determine the lightness of an image, 

Following a different approach, Ullman(1979a) studied the perception of apparent motion generated 
by two successive frames consisting of isolated dots of equal intensity moving independently of each other. 
Without constraint, all possible pairings, or "correspondences", of dots in the first frame with dots in the 
second are equally likely. Ullman defined the "most likely" correspondence to be the one that minimized the 
sum 

!<*<n 
!</<m 

where n is the number of dots in the first frame, m is the number of dots in the second frame, and sy is one if 
the ith dot of the first frame P { is paired with the jth dot of the second frame Q 3i else zero. The weight q^ is 
the "cost" of pairing P, with Q J9 and might, for example, be related to the image distance between the paired 
points. The problem of finding the minimal correspondence is considered in terms of integer programming. If 
correspondences are assumed to be covering mappings, the following linear constraints apply to the xtf. 



Vi, 1 < % < n ^2 Xij > 1, 

l<j<m 

and 



Vj,l<i<m ]H *ij>l. 



Ki<n 
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Ullman restricted the set of Qj that can be paired with Pi to those whose positions were close to P t . Following 
Arrow, Hurwicz, and Uzawa(1958), he set up the iterative scheme 



l<t<n 
l<i<m 

The approach can be extended to mappings that are not one-one, as well as to continous motion. A major 
problem with the approach is guaranteeing the convergence of the algorithm. This is determined largely by the 
properties of the costs g y , but this was not investigated, aside from a comment on the empirical determination 
of the qij (see also Ullman 1979b). 

One limitation of Ullman's approach is that it is restricted to minimizing a known linear objective func- 
tion that is subject to linear constraints. The method can be extended to constrained nonlinear programming 
in which the goal is to minimize a known function f(x) subject to a set of equality and inequality constraints 
of the form &(s) < 0. In general, however, criteria based on other than intuition need to be found for 
selecting the function / to be minimized. To do this, one can apply the calculus of variations (see for example 
Courant and Hilbert 1953, chapter IV). The familiar differential shows how to find a real valued parameter 
that minimizes some function. The calculus of variation extends the differential calculus by showing how one 
can determine a junction f \ which is subject to a given set of boundary conditions, and minimizes the integral 

*(/) = J J ' F[*,vJ,fxtfyJzx>fxy,f w )dxdy (17) 

over a given region of integration G\ The function F is called a "performance index" and generalizes the 
notion of cost function associated with linear and nonlinear programming. In the next section we shall con- 
sider the choice of a performance index for interpolating smooth surfaces from one-dimensional boundary 

conditions. 

t For simplicity of presentation, we restrict attention to functions / of one or two variables x, y, 



J0*\_ 
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Associated with a variational problem of the form (17) is the Euler equation, which provides a necessary, 
though by no means sufficient, condition which a function / must satisfy if it is to minimize the variational 
integral SF(/). For the particular variational problem given in equation (17), the Euler equation is 



Ff ~ dx F} * ~ ¥y F ^ + M Ff ** + W^ + V F/vv = °' (18) 

In the case that there is only a single dependent variable x 9 the partial derivatives are total and the Euler 
equation becomes 



F >-i F '-+£? F >~=< > - < 19 > 



There are two important considerations associated with the use of the calculus of variations. First, unlike 
the differential calculus, the existence of an extremum f of the integral given in equation (17) cannot be taken 
for granted. Courant and Hilbert(1953, p. 173) note that "a characteristic difficulty of the calculus of variations 
is that problems which can be meaningfully formulated may not have solutions". Conditions for the existence 
of a minimum have recently been proposed by Grimson(1981) and will be discussed in the next section. 

Second, associated with any variational problem is a set of natural boundary conditions which imposes a 
necessary condition on any feasible solution to the Euler equation at the boundary. Courant and Hilbert(1953, 
p. 211) note that "in general, we can, by adding boundary terms or boundary integrals essentially modify 
the natural boundary conditions without altering the Euler equations". Determining the "most conservative" 
solution means finding a performance index that guarantees the existence of an extremum function /* and 
provides the tightest set of natural boundary conditions that are consistent with the given data. 

The calculus of variations has recently been applied by a number of authors to interpolate plane and 
space curves and surfaces. We review the applications in that order. First, Horn(1981) has recently determined 
the curve which passes through two specified points with specified orientation while minimizing 
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K 2 ds, (20) 



where k is the curvature and a is the arc length. This is the true shape of a spline used in "lofting" (Faux and 
Pratt 1979,p. 228). In a thin beam, curvature is proportional to the bending moment. The total elastic energy 
stored in a thin beam is therefore proportional to the integral of the square of the curvature. Since the shape 
taken on by a thin beam is the one which minimizes the internal strain energy, the curve that solves equation 
(20) is called the "curve of least energy". The variational problem is to minimize 



/ 



rd$. 



(i + /2)» 

This has the form of equation (17). Horn(1981, page 19) shows that the Euler equation is 



— CK = \/ cos i>* 

where t/> is the angle between the tangent to the curve and the axis of symmetry. The solution to this 
differential equation is an incomplete elliptic integral of the first kind. Brady, Grimson, and Langridge(1980) 
consider a "small angle" approximation to the curve of least energy, in which first derivatives can be ignored. 
The performance index that they use is f* x , for reasons that will become evident in the next section. They find 
that in that case the solution is a cubic. Horn(1981,page 2) notes that the fact that a curve has near minimum 
energy does not mean that it lies close to the curve of minimum energy. Note that the existence of the curve 
of least energy is guaranteed as Horn has derived an analytical formula for it. Approximations to it, such as 
Brady, Grimson, and Langridge's are similarly guaranteed to exist. 

Barrow and Tenenbaum(1981) investigate the problem of interpreting a line as the image of a space curve 
that is an occluding boundary. They observe that the problem has two parts: (i) determining the tangent 
vector t at each point on the space curve, and (ii) determining the surface normal at each point, given that it is 
constrained to be orthogonal to the tangent. 



^•*K_ 
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They suggest minimizing a performance index F that is a function of the curvature k and the torsion r 
(possibly together with their derivatives), and expresses a suitable notion of "smoothness". They first consider 
uniformity of curvature as a measure of smoothness, that is F = gj = k u where s measures distance along 
the space curve. They reject this measure on the grounds that/ci can be made arbitrarily small by "stretching 
out the space curve so that it approaches a twisting straight line". To overcome this difficulty, they propose 
that the space curve should also be "as planar as possible or, more precisely, that the integral of its torsion 
should be minimized". 

Barrow and Tenenbaum finally suggest finding the space curve that projects to the given image line and 
minimizes the performance index [^^] 2 , where b is the binomial. They report that an algorithm based on 
their analysis produced the "correct 3-D interpretations for simple and closed curves, such as an ellipse, which 
was interpreted as a circle". However, they note that the rate of convergence was slow and dependent on the 
initial data. No consideration is given to the Euler equations, to the existence of an extremum given a line 
drawing {x(s) } y(s)}, or to the natural boundary conditions associated with the performance index [~j~] 2 . 
Empirical evidence that the method works on a number of simple test cases is encouraging; but there is no 
analysis of the scope of the method. 

In the same paper, Barrow and Tenenbaum(1981) consider the interpolation of a smooth surface from 
depth and local surface orientation values at all points along the surface boundary. Their approach is to 
"seek a technique that yields exact reconstructions for the special symmetric cases of spherical and cylindrical 
surfaces, as well as intuitively reasonable reconstructions for other smooth surfaces." (Barrow and Tenenbaum 
1981). They observe that if n is the surface normal of a cylinder, then the x and y components of the normal 
n x and n y are linear functions of x and y, so long as the axis of the cylinder lies in the x — y plane. This 
observation forms the basis of an algorithm to estimate the surface normal by least squares fitting of the 
parameters of the partial derivatives of the normal. As before, no analysis is given of the Euler equation, the 
natural boundary conditions, nor the convergence of their algorithm for different types of surface. 
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5. A performance index for surface interpolation. 

In the review of the application of the calculus of variations to visual perception in the previous section 
we drew attention to three important considerations. First, the Euler equations provide a necessary condition 
on possible extremal functions. Second, the existence of an extremum cannot be taken for granted, even when 
the minimization problem seems plausible on some grounds. Third, the natural boundary conditions impose 
a necessary condition on any feasible solution to the Euler equation at the boundary. The most thorough 
analysis of the second of these problems in Computer Vision, framed in the context of surface interpolation, is 
due to Grimson(1981), who proves the following theorem. 

Theorem (Grimson, see Rudin(1973)) Suppose there exists a complete semi-norm F on a space of func- 
tions 9, and that F satisfies the parallelogram law. Then, every non-empty closed convex set 6C? contains a 
unique element/* of minimal normF(/*), up to possibly an element of the null space ofF. 

A semi-norm F is a function Vh»+ from a vector space V to the positive real numbers that satisfies 

F(t> + u>)<F(v)+F(u>) 
F(av) = \a\F{v). 

Informally, a semi-norm is a generalization of the Euclidean metric, and provides a measure of a vector. The 
second condition generalizes the triangle inequality, for example. The null space of die semi-norm F consists 
of all those vectors vo that map to zero. Since 

F(v + vo)=F(v), 

any element of the null space can be added to a vector of minimal norm to yield another vector of minimal 
norm. Hence the qualifying phrase "unique ... up to possibly an element of the null space of F". The 
parallelogram law states that 

[F(v + w)} 2 + [F(v - w)} 2 = 2[F(v)] 2 + 2[F(w)} 2 , 



^*\ 
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for all vectors v, w. Finally, the semi-norm is complete if all Cauchy sequences converge. As is well known, 
the elements of vector spaces can be functions. This enables Grimson to prove the following Corollary, that 
guarantees the existence of an extremum function in calculus of variations "most conservative" interpolation 
problems. 

Corollary (Grimson 1981). Let the set of known points be {(xi, y*, Zi) | 1 < t < n}. Let # be a vector 
space of possible functions on 5R 2 and let 6 be the subset of ^ that interpolates the known data. That is, for all 
functions /e8, f(xi t yi) = Z{. Let F be a complete semi-norm on 8 that satisfies the parallelogram law. Then 
there exists a ynique (up to the null space of F) function /* that interpolates the data and has minimal norm. 
In particular, if F is a performance index then there is a function f that minimizes the integral 



m = fi 



In short, if the conditions of the Corollary are fulfilled, the existence of a "most conservative" surface that 
meets the boundary conditions is guaranteed. As we shall see, the condition of being a semi-norm is the most 
restrictive required of the performance index. The conditions are sufficient to guarantee the existence of a 
minimum, but they are not necessary. For example, /c 2 is not a seminorm*; nevertheless Horn*s(1981) analysis 
shows that there is a unique minimum. It is far from clear whether Barrow and Tenenbaum's(1981) analyses of 
curve and surface interpolation have a guaranteed minimum in all cases. 

Grimson notes that several intuitively plausible performance indices are not semi-norms. For example, 
the two most popular measures of curvature are not. Suppose that k\ and /C2 are the principal curvatures of 
a surface(Faux and Pratt 1979, p. Ill), then the Gaussian curvature k 9 is the product k^ and the mean 
curvature K m is the sum /ci + «2. For a surface f(x t y), 



, . fxxfyy J X y 



(l+fl + Py) 2 



/~\ 



V/hich is why Brady, Grimson, and I^ngridge(1980) used the small angle approximation f% x . 



28 



Since the curvatures can be negative, while a semi-norm is required to be positive, it is necessary to 
investigate 



\J K]dxdy\ 



Grimson(1981) observes that/c^(a/) =^ |'a|«J(/) because of the denominator. If f x and f y are small, the 
denominator is approximately equal to one, and the numerator is a seminorm. Note that it is 

fxxfyy—fxy (21) 

Grimson shows that the mean curvature K m is also not a semi-norm for exactly the same reason. The 
analogous small angle approximation is 

(/,x + /yy) 2 -(A/) 2 , 

the square Laplacian, which is a semi-norm. We find it convenient to denote the square Laplacian by JFJ. 
Grimson(1981) chooses the quadratic variation 

J xx i **! %y I Jyyf 

on the grounds that its null space, consisting of all linear functions, is smaller than the null space of the square 
Laplacian. If we denote -.the quadratic variation by F qi we see that the approximation to the Gaussian curvature 
given in equation (21) is * '7 • 

How shall we choose a performance index for surface interpolation, given that it has to satisfy the condi- 
tions of the Corollary? We have exhibited three candidates, are there more? Notice first that each of the 
semi-norms given above are quadratic forms in f XXi f xy , and f yy . It is easy to show that any quadratic form 
satisfies the semi-norm and parallelogram conditions, and so there, is an infinite set of plausible semi-norms to 



. /■"■^ 
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use to find the "most conservative" interpolated surface. We need an extra condition, and the one we choose 
is rotational symmetry, since we suppose that surface interpolation is an isotropic process. Proposition 6 of 
section 3 shows that the rotationally symmetric quadratic forms in f xx , f xy9 and f yy form a vector space that 
has the square Laplacian and the quadratic variation as a basis. The choice of which performance index to use 
is thus effectively reduced to the square Laplacian, the quadratic variation, and linear combinations of them. 
How shall we choose between those two? In the light of our earlier discussion, two criteria suggest themselves: 
the Euler equations and the natural boundary conditions. 

Proposition 7. All rotationally symmetric quadratic forms lead to an identical Euler equation 

A 2 (/) = 0. 

Proof. We exploit the fact that the square Laplacian and the quadratic variation are a basis of the 
rotationally symmetric quadratic forms. 

a.Square Laplacian: The performance index is 

F, = (/^ + / yv ) 2 . 

By equation (18) the Euler equation is 



j^{2(/** + fyy)} + j^{2(/** + fyy)} - 0, 



that is 



(A/) 2 = 0, 



as required. 

b. Quadratic variation: The Euler equation is 
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"fxzzx "T *fxyxy "T ^fyyyy ' — "» 



that is 



(A/) 2 = 0, 

provided that / is continuous of fourth order. 

c.Linear combinations ofF\ andF q .\ Linear combinations clearly give rise to the identical Euler equation^ 
The gist of Proposition 7 is that there is no difference between F q and Fi in the interior away from the 

boundary conditions. We can see the result of Proposition 7 in an alternative interesting way. Recall that 

Fi l ~— Fq = 2{fxxfyy ™ f X y)t 

is the semi-norm approximation to the Gaussian curvature (equation 21). The latter expression is an instance 
of a divergence expression, and Courant and Hilbert(1953, p. 196) note "If the difference between the in- 
tegrands of two variational problems is a divergence expression, then the Euler equations and therefore the 
families of extremals are identical for the two variational problems/ 1 

Since F q and JF] have identical Euler equations, we analyze their natural boundary conditions in order 
to choose between them. We could approach this problem directly; but a more revealing route is available. 
Courant and Hilbert(1953, p250) consider the statics of a thin plate. In particular they determine the shape it 
assumes for a given force p{s) along its boundary T and bending moments m(s) normal to its boundary. 

Courant and Hilbert note that the energy stored in the plate is the integral of a quadratic form in the 
principal curvatures /ci and k 2 of the surface, a result which can be derived from noting that the elastic energy 
stored in a thin strip (corresponding to any normal section) is proportional to the square curvature. It follows 
that the stored energy is locally 



31 



^\ 



/**\ 



g 1 = Q (Acf + k\) + 2/3/ci/ci 

== o(/ci + «2) 2 + 2(/? — Q)/CiKi 

= CWC m + 2(/? — Q>C ff 

= /?F< + (a-/3)F g 
= a(MF, + (l-M)P ff ), 

where ^ = -. It follows that the energy stored in the thin plate is a convex linear combination of the 
square Laplacian and the quadratic variation, which formally establishes its connection to the visual percep- 
tual problem studied here. Observe that setting the weight \i — 1 gives the square Laplacian, while setting it 
equal to zero gives the quadratic variation. Note also that this expression for the stored energy makes use of 
the small angle approximation to the curvature used in equation 21. 

A second source of stored energy derives from the boundary conditions that are represented as a function 
p(s) along the boundary T of the plate and a bending moment m(s) applied normal to the plate. Courant and 
Hilbert(1953, p. 251) show that the natural boundary conditions associated with the plate are 

-A/ + (1 - fj){f xx x 2 8 + 2f xy x 8 y a + f yy yl) = p{s) 
A/ + (1 — l*)^-{fxxX g x n -f f xy {x a y n -f x n y s ) + fy y y a y^ = m(s), 



dn 



that is 



'da y 



-A/ + (1 - n)({*»Vs)H[x 3 y a ) T ) = p{ 8 ) 
JU/ + (1 - /^([z„y n ] H [x a y a ) T ) = m(s), 



where H is the Hessian matrix 



Jxx Jxy 

Jxy fyy 

Gladwell and Wait(1979) quote version of this result due to Agmon(1965), that the biharmonic operator, 
which we showed was the natural boundary condition for the surface interpolation problem, has Dirichlet 
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forms that are linear combinations of the square Laplacian and the quadratic variation. As an example of the 
constraint, consider a straight line contour aligned with the z-axis. Then [x 8 y a ] — [1 0] and [x n y n ] = [01]. 
The natural boundary conditions reduce to 

fyyy + (2 ~ l*)fyx% = m(a). 

The constraint is tightest when \i is not equal to one. A similar result can be obtained for a straight line 
contour inclined at an angle a to the s-axis. The first of the natural boundary conditions is 

/^(sin 2 a + \x cos 2 a) -f / y2/ (cos 2 a + /i sin 2 a) + (1 — /*) sin 2a/ xy . 

If fi = 1, there is no constraint from the cross derivative. If// is not equal to 1, at most one of the 
terms can be zero. We conclude that interpolation problems in which the small angle approximations used 
throughout our analysis hold it is preferable to choose \x not equal to one, that is to say to not use the square 
Laplacian as a performance index. The quadratic variation is an obvious choice, but so are linear combinations 
of the square Laplacian and the quadratic variation for which /i is not equal to one. Grimson(1981) chooses 
the quadratic variation since its null space is smaller than that of the square Laplacian. This is a precise way 
of saying that it imposes a tighter constraint. For example, the function f(x, y) = xy is in die null space of 
the square Laplacian but not in the null space of the quadratic variation. Since the quadratic variation has 
the smallest null space among the linear combinations of the square Laplacian and quadratic variation, this 
is an additional reason for choosing it. We would further expect that any differences between the quadratic 
variation and the square Laplacian would show up near the given boundary data but not in die interior, far 
removed from the boundary. This is what Grimson(1981) finds in a set of examples that compare surfaces 
interpolated using the quadratic variation and the square Laplacian. 
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